Graphene, a monolayer carbon sheet, possesses many excellent properties, such as high electric conductivity,^[1–4] high thermal conductivity,^[5–7] and high stability,^[8,9] which results in extensive applications in various fields. By stacking graphene, two- or few-layer graphene can be formed. Beside the conventional AA and AB stacking graphene, different stacking sequences can result in various Moire patterns between layers.^[10–12] In experiment, these twisted few-layer graphenes (T-FLGs) have been synthesized by using the chemical vapor deposition (CVD), mechanical exfoliation or growth on the carbon terminated SiC surfaces.^[13–15] Although the interlayer interactions are weak Van der Waals force, the T-FLG shows completely different electronic properties from those of AA and AB stacking graphene structures, and their properties are strongly dependent on the rotation angle θ. For example, Bstristzer and MacDonald studied the electronic structure of a twisted bilayer graphene, and found that in the continuum Dirac model, the Moire pattern periodicity leads to Moire Bloch bands;^[16] Li et al. demonstrated that a rotation between graphene layers can generate Van Hove singularities, which can be brought arbitrarily close to the Fermi energy by varying the angle of rotation;^[17] Dos Santos et al. also found that a significant reduction of Fermi velocity occurs in twisted graphene. Moreover, unlike the case of the AB stacked bilayer, a potential difference between layers does not open a gap in the spectrum.^[18]

Recently, phonon properties of T-FLG have received attention gradually. Several experimental Raman studies found a series of phonon peaks in different energy regions, which is absent in the few-layer graphene without twist.^[19,20] Cocemasov et al. found that the twisted bilayer graphene leads to the emergence of new phonon branches, termed entangled phonons.^[21] Nika et al. found that the phonon specific heat has an intriguing dependence on the rotation angle θ in twisted bilayer graphene, which is particularly pronounced at low temperature.^[22] These studies indicate that the twisted graphene structures have unique thermal properties. To the best of our knowledge, however, thermal transports in these twist graphenes have been reported little. Therefore, it is natural to ask what is the effect of twist on the thermal transport in few-layer graphene.

In this paper, we investigate the thermal transport properties of T-FLG, by using molecular dynamics (MD) simulations. Thermal transports in the parallel and perpendicular directions are considered, respectively. Our results show that perpendicular and parallel thermal conductivities κ_⊥ and κ_|| are both strongly dependent on rotation angle θ, and θ = 30° is a transition angle. Underlying mechanisms for the unique phenomenon are explained by out-of-plane phonon density of states (PDOS) and interlayer interaction energy.

Figures 1(a) and 1(b) show the simplest models for perpendicular and parallel transport in T-FLGs, respectively. Figure 1(a) is a three-layer structure, where hot and cold heat baths are placed on the top and bottom layers, respectively. Figure 1(b) shows a bilayer graphene where the hot and cold baths are located on the left and right terminals, respectively. The stacking forms of the T-FLGs can be tuned by rotating one of the graphene layers. Figure 1(c) shows an initial twisted bilayer graphene, which is in fact an AB stacking structure with a rotation angle θ = 0°. The red dot in the red circle is denoted as a rotation center. When the top layer of graphene rotates clockwise, various twisted structures with different Moire patterns can be formed to be dependent on θ. Figure 1(d) presents twisted structures with θ = 0°, 15°, 30°, 45°, and 60°, respectively. One can note that the rotating period is 60°, in other words, the structure with 60° + θ is the same as the structure with θ, while θ = 60° corresponds to an AA stacking bilayer graphene. In the θ = 30° structure, the zigzag edges of the top layer align with the armchair edges of the bottom layer. While the patterns of the structures with θ = 15° and 45° are more complicated. For the perpendicular transport, the sizes of all graphene layers are fixed at width W = 6.31 nm and length L = 6.24 nm. For the parallel transport, the width and length of graphene layer are W = 5.02 nm, and L = 13.01 nm. To study the effect of twist on the thermal transport in T-FLG, twisted structures with more than 2 and 3 layers are considered. In these structures, the stacking form between even layers is AA as well as that between odd layers, while the stacking between even and odd layers can be twisted by a rotation angle θ as shown in Figs. 1(e) and 1(f).

Fig. 1.

Figure Option ViewDownloadNew Window

Fig. 1. (color online) Schematic views of T-FLGs. (a) Simplest model for perpendicular transport in twisted three-layer graphene. The red and blue regions represent hot and cold heat baths, respectively. (b) Simplest model for parallel transport in twisted bilayer graphene. (c) An initial twisted bilayer graphene with a size of W × L, which is an AB stacking bilayer graphene with rotation angle θ = 0°. The red dot represents the rotation center. Other twisted graphenes can be formed by clockwise rotating the top layer. (d) Twisted bilayer graphene structures with θ = 0°, 15°, 30°, 45°, and 60°, respectively. (e) Top view of a T-FLG with a rotation angle θ. Stacking form between even (or odd) layers is AA, while θ is the rotation angle between even and odd layers. (f) The side view of a T-FLG with a rotation angle θ.

Thermal transport of T-FLG is calculated by MD method. All MD simulations here are performed by using the LAMMPS package.^[23,24] An optimized Tersoff potential^[25,26] is used to describe the interactions between intra-layer carbon atoms. The 12-6 Lennard–Jones (LJ) potential is adopted to describe the interlayer van der Waals (VDWs) interactions:^[27] 1 where r represents the distance between two atoms, ε is the energy that reflects their interaction strength, and σ denotes the zero-across distance of the potential. The parameters of the LJ potential are cited from Refs. [28] and [29]. The cutoff distance of VDWs force is chosen to be 10 Å. We first relax all the twisted structures under the canonical ensemble (NVT) conditions for 50 picoseconds (ps), in time steps of 0.5 femtoseconds (fs). The heat baths are Nosé–Hoover thermostats, and the temperatures of hot and cold baths are T₀ + ΔT and T₀ − ΔT, respectively, to achieve a temperature gradient. T₀ is set to be environment temperature (300 K) and the temperature drop is set to be ΔT = 30 K. In order to avoid the graphene moving in the simulations, we use fixed boundary conditions at the two ends of the graphene.^[30,31] From Fourier’s law, we define the thermal conductivity κ as^[32] 2 where ∇T is the temperature gradient, J is the heat flux from the heat source to sink, and S is the cross-sectional area. To reach the non-equilibrium steady states, the system is relaxed for 50 ps under NVE ensemble. After that, 60 ps is used to calculate the thermal conductivity. The previous studies showed that using optimized Tersoff potential and LJ potential parameters can well describe the phonon properties and thermal transport behaviors in few-layer graphene, indicating that our calculated results are accurate and reliable.^[33–35]

We first consider perpendicular thermal transport in the T-FLGs of Fig. 1(c), where heat flux flows across three graphene layers along the perpendicular direction. The top and bottom layers are AA stacking, while the middle layer is rotated to tune the heat flux. Figure 2(a) presents thermal conductivity κ_⊥ as a function of rotation angle θ. We find from Fig. 2(a) that κ_⊥ exhibits a non-monotonic dependence on θ. When θ is changed from 0° to 60°, κ_⊥ first decreases gradually and then increases sharply. The lowest κ_⊥ occurs at θ = 30°, which is a transition point. The thermal conductivity at θ = 30° is approximately 3.5 × 10⁻³ W/mK, while those at θ = 0° and 60° are 4.3 × 10⁻³ and 8.0 × 10⁻³ W/mK. The latter are 1.3 and 2.3 times the former, respectively. All these values are small because of weak interlayer interactions between graphene layers. However, one can find that the perpendicular transport is sensitive to the stacking form and rotation angle.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (color online) (a) Perpendicular thermal conductivity κ⊥ of T-FLG in Fig. 1(a) as a function of rotation angle θ. (b) Comparison among out-of-plane PDOSs for three T-FLGs with θ = 0°, 30°, and 60°. Inset shows the comparison among in-plane PDOSs for three T-FLGs with θ = 0°, 30°, and 60°. (c) Comparison among interlayer interaction energies of the three T-FLGs.

In order to explain the effect of twist on the thermal transport, we compare the out-of-plane PDOSs of the three-layer graphene structures with θ = 0°, 30°, and 60° as shown in Fig. 2(b). One can find that the PDOSs of the three structures exhibit different spectra, especially at the low phonon frequencies (< 10 THz). The PDOSs of θ = 60° and 30° present the highest and lowest values, respectively, while that of θ = 0° is in the middle. It is known that the out-of-plane phonon mode ZA of graphene is in a frequency range of 0–10 THz, and the perpendicular transport is strongly dependent on the ZA mode.^[36–38] Therefore, the higher PDOS induces higher κ_⊥. In order to further explore the relation between interlayer interactions and rotation angle θ, the interaction energies of the three structures are compared in Fig. 2(c). It is found that the interaction energy of θ = 60° is largest, while that of θ = 30° is the smallest. In the case of θ = 60°, all layers are of AA stacking, and atoms between layers are directly stacked on each other. The shortest distance between atoms induces strong interlayer interaction. In the case of θ = 0°, the stacking form between two nearest layers is AB, where only half of the atoms are directly stacked. In the case of θ = 30°, the number of atoms having direct shortest-distance interactions is the smallest. Therefore, the interaction energy of θ = 60° is the largest, while that of θ = 30° is the smallest. The weaker interaction leads to weak coupling of out-of-plane phonon modes, and thus the thermal conductivity κ_⊥ is lower.

In Fig. 3, the effect of layer number on κ_⊥ of T-FLGs is calculated. The black, red, and blue lines in Fig. 3(a) present κ_⊥ values of T-FLGs with θ = 0°, 30°, and 60° each as a function of the number of layers, respectively. It is found that all the thermal conductivities increase with layer number increasing gradually, and then converge to steady values after 13 layers. This indicates that thermal transport in the few-layer structures approaches the transport in the bulk gradually after 13 layers. Relatively, the thermal conductivity of θ = 60° grows fast as the layer increases. The steady value is approximately 35.0 × 10⁻³ W/mK, which is 3 times more than that of θ = 30° structure. The difference between the structures of θ = 0° and 30° seems small. To clearly compare the difference between θ = 0° and 30° structures, the values of normalized heat flux J = J_N/J₀ of the two structures each as a function of layer number are shown in Fig. 3(b). Here J_N is the heat flux of twisted N-layer structure, while J₀ is the heat flux of three-layer graphene with θ = 0°. The heat fluxes are calculated by setting heat baths separately on top and bottom layers, and temperature difference ΔT = 30 K is fixed. With the increase of the layer number, the heat fluxes of θ = 0° and 30° both first drop and then converge gradually. One can note that the trend of heat fluxes in Fig. 3(b) is inverse to that of thermal conductivities in Fig. 3(a). This is because the decrease of temperature difference between layers is faster than the increase of thermal conductivities when the layer number is increased. One can clearly find that the heat flux in θ = 30° structure drops more rapidly than that in θ = 0° structure. From these results, one can conclude that the T-FLGs inherit the thermal transport characteristics of the three layer structures: the θ = 60° structure has the highest thermal transport ability while θ = 30° structure has the weakest ability. Because all the T-FLGs possess small thermal conductivities, i.e., high thermal resistances, these materials have great potential applications in thermal insulation.^[39–41] One can tune the insulation effect by rotating the angle between layers. Obviously, the structure with θ = 30° exhibits the best insulation effect.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. (color online) (a) Plots of perpendicular thermal conductivity κ⊥ of T-FLGs versus the number of layers. (b) The comparison between the values of normalized heat flux J = JN/J0 for T-FLGs with θ = 0° (black dash-dotted line) and 30° (red dash-dotted line), where JN is the heat flux of twisted N-layer structure and J0 is the heat flux of three-layer graphene with θ = 0°.

Next, we consider parallel thermal transport in T-FLGs where heat flux flows in the structure along the parallel direction. In Fig. 4(a), we calculate the values of parallel thermal conductivity κ_|| of twisted two-layer graphene in Fig. 1(b) with different values of rotation angle θ. One can find that, the relation between thermal conductivity and rotation angle is completely different from that in the case of perpendicular transport. Here, κ_|| of θ = 30° has the highest value, while those of θ = 0° and 60° are smaller. The parallel transport is related to two processes: one is the transport in the graphene layers, the other is the phonon scattering between layers because of the interlayer coupling. Because the layer coupling is relatively weak, the direct effect of interlayer interaction on the parallel transport is very limited. This can be illustrated by comparing the in-plane PDOSs of three structures, which is shown in the inset of Fig. 2(b). Therefore, the parallel transport is dominated by the edges of graphene layers. For the θ = 30° structure, along the transport direction, the two graphene layers have zigzag and armchair edges, respectively. For the θ = 0° and 60° structures, both layers have the armchair edges here. The thermal conductivity of zigzag nanoribbon is larger than that of armchair nanoribbon because of anisotropic transport in graphene.^[42,43] Therefore, the κ_|| value of θ = 30° is larger than those of θ = 0° and 60°. To describe the effect of interlayer scattering on the thermal transport, in Fig. 4(a) the κ_|| values of the structures without interlayer coupling are also presented for comparison, and the difference Δκ_|| between the structures with and without coupling is shown in Fig. 4(b). It is seen from Fig. 4(b) that the difference of thermal conductivity is sensitive to the rotate angle, and exhibits a somewhat similar trend to that of the κ_⊥ in Fig. 2(a). This indicates that the interlayer scattering in the θ = 30° structure is the smallest because of the weak interlayer interaction mentioned above. Therefore, the strength of interlayer coupling has an inverse effect on the perpendicular transport compared with that on the parallel transport. It is conducible to the perpendicular transport while suppressing the parallel transport.

Fig. 4.

Figure Option ViewDownloadNew Window

Fig. 4. (color online) (a) Parallel thermal conductivity κ|| of twisted bilayer graphene in Fig. 1(b) as a function of rotation angle θ, where the black and red lines represent the cases with and without interlayer coupling. (b) Difference between thermal conductivities Δκ|| (= κ|| without coupling −κ|| with coupling) varying with rotation angle θ. (c) Variations of κ|| of T-FLGs with θ = 0°, 30°, and 60° with the number of layers. (d) The values of Δκ|| of T-FLGs with the number of layers.

In Fig. 4(c), we calculate κ_|| values of three T-FLGs (θ = 0°, 30°, and 60°) varying with layer number. All κ_|| values decrease rather than increase with layer number increasing, which is inverse to the trend in Fig. 3(a). This is because more graphene layers strengthen phonon scattering and thus weaken the thermal transport along the parallel direction. With the increase of layer number, the interlayer scattering becomes steady gradually and then the thermal conductivities converge. To compare the interlayer scatterings in the three structures, we calculate the difference Δκ_|| between the cases with and without interlayer coupling in these structures as shown in Fig. 4(d). The Δκ_|| of θ = 60° increases rapidly with the increase of layer number, while that of θ = 30° increases slowly. This further proves that the θ = 30° structure has the weakest interlayer interaction. Li et al. investigated experimentally the thermal conductivity of suspended twisted bilayer graphene by using an optothermal Raman technique, and they found that the thermal conductivity of twisted bilayer graphene is lower than those of monolayer graphene and AB stacked bilayer graphene in the entire temperature range examined (∼ 300–700 K).^[44] These are consistent with our results.

In this work, we study perpendicular and parallel thermal transports in the T-FLGs, by using MD simulations. The effect of twist on the thermal transport is discussed by comparing the thermal conductivities κ_⊥ and κ_|| in the T-FLGs with different values of rotation angle θ. When θ is changed from 0° to 60°, κ_⊥ first decreases and then increases, where θ = 30° is the transition point. This indicates that T-FLG with θ = 30° has the strongest interlayer coupling, which is proven by analyzing PDOS and interlayer interaction energy. For the parallel transport, the relation between thermal conductivities κ_|| and θ is complicated, because intra-layer thermal transport is more sensitive to the edge of the layer than their stacking form. However, the dependence of interlayer scattering on θ is similar to that on the κ_||, which reveals the inverse effect of interlayer coupling on perpendicular and parallel transports. The increase of layer number further increases the twist effect in the T-FLG. These results provide an insight into understanding thermal transports in two-dimensional VDWs materials.